﻿/******************************************************************************
 * 
 * Announce: CSharpKit, Basic algorithms, components and definitions.
 *           Copyright (C) ShenYongchen.
 *           All rights reserved.
 *   Author: 申永辰.郑州 (shenyczz@163.com)
 *  WebSite: http://github.com/shenyczz/CSharpKit
 *
 * THIS CODE IS LICENSED UNDER THE MIT LICENSE (MIT).
 * THIS CODE IS PROVIDED *AS IS* WITHOUT WARRANTY OF 
 * ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING ANY
 * IMPLIED WARRANTIES OF FITNESS FOR A PARTICULAR
 * PURPOSE, MERCHANTABILITY, OR NON-INFRINGEMENT.
 * 
******************************************************************************/

using System;
using CSharpKit.Numerics.Theories;

namespace CSharpKit.Numerics.LinearAlgebra.Factorization
{
    /// <summary>
    /// LU - LU分解
    /// </summary>
    /// <remarks>
    /// For a matrix A, the LU factorization is a pair of lower triangular matrix L and
    /// upper triangular matrix U so that A = L*U.<para/>
    /// In the Math.Net implementation we also store a set of pivot elements for increased
    /// numerical stability. The pivot elements encode a permutation matrix P such that P*A = L*U.<para/>
    /// </remarks>
    /// <typeparam name="T"></typeparam>
    public abstract class LU<T> : ISolver<T>
        where T : struct, IFormattable, IEquatable<T>
    {
        #region Constructors

        protected LU(Matrix<T> factors, int[] pivots)
        {
            Factors = factors;
            Pivots = pivots;

            _lazyL = new Lazy<Matrix<T>>(ComputeL);
            _lazyU = new Lazy<Matrix<T>>(Factors.UpperTriangle);
            _lazyP = new Lazy<Permutation>(() => Permutation.FromInversions(Pivots));
        }

        #endregion


        static readonly T One = Numeric<T>.One;


        #region Fields

        readonly Lazy<Matrix<T>> _lazyL;
        readonly Lazy<Matrix<T>> _lazyU;
        readonly Lazy<Permutation> _lazyP;

        protected readonly Matrix<T> Factors;   // L
        protected readonly int[] Pivots;        // P => PA=LU

        #endregion


        #region Properties

        public Matrix<T> L => _lazyL.Value;
        public Matrix<T> U => _lazyU.Value;
        public Permutation P=> _lazyP.Value;

        /// <summary>
        /// 行列式
        /// </summary>
        public T Determinant => _GetDeterminant();

        #endregion


        /// <summary>
        /// 
        /// </summary>
        /// <returns></returns>
        public abstract Matrix<T> Inverse();



        #region ISolver<T>

        public Matrix<T> Solve(Matrix<T> input)
        {
            var x = Matrix<T>.Builder.SameAs(input, input.RowCount, input.ColumnCount, fullyMutable: true);
            Solve(input, x);
            return x;
        }
        public abstract void Solve(Matrix<T> input, Matrix<T> result);

        public Vector<T> Solve(Vector<T> input)
        {
            var x = Vector<T>.Builder.SameAs(input, input.Count);
            Solve(input, x);
            return x;
        }

        public abstract void Solve(Vector<T> input, Vector<T> result);

        #endregion



        #region private functions

        /// <summary>
        /// 计算下三角矩阵
        /// </summary>
        /// <returns></returns>
        private Matrix<T> ComputeL()
        {
            var result = Factors.LowerTriangle();

            for (var i = 0; i < result.RowCount; i++)
            {
                result.At(i, i, One);
            }
            return result;
        }

        private T _GetDeterminant()
        {
            var det = Numeric<T>.One;

            for (var j = 0; j < Factors.RowCount; j++)
            {
                if (Pivots[j] == j)
                {
                    //det *= Factors.At(j, j);
                    det = NumOperators.Multiply(det, Factors.At(j, j));
                }
                else
                {
                    //det *= -Factors.At(j, j);
                    det = NumOperators.Multiply(det, NumOperators.Negate(Factors.At(j, j)));
                }
            }

            return det;
        }

        #endregion


        //}}@@@
    }


}

